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Abstract 

•i-H 

X. 

H , We present Monte Carlo simulation results for the three dimensional Thirring 

model on moderate sized lattices using a hybrid molecular dynamics algorithm which 
permits an odd or non-integer number Nf of fermion flavors. We flnd a continuous 
chiral symmetry breaking transition for Nf ^ 3 with critical exponents consistent 
with expectations from previous studies. For Nf = 5 the order of the transition is 
difficult to determine on the lattice sizes explored. We present a phase diagram for 
the model in the {l/g'^,Nf) plane and contrast our flndings with expectations based 
on approximate solutions of the continuum Schwinger-Dyson equations. 

PACS: ll.lO.Kk, 11.30.Rd, 11.15.Ha 

Keywords: four-fermi, Monte Carlo simulation, dynamical fermions, chiral symmetry 
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1 Introduction and Algorithm 

Recent numerical investigations [|l], @, @] have revealed an interesting phase structure 
for the three dimensional Thirring model as a function of coupling g"^ and number of 
fermion species Nf. For sufficiently small Nf < Nfc, there is a continuous transition 
at g"^ = g'^{Nf) between a weak coupling phase in which chiral symmetry is realised in 
the limit m ^ 0, and a strong coupling one in which the symmetry is spontaneously 
broken. The critical indices characterising the transition are distinct for Nf = 2 and 
Nf = 4, suggesting that the continuum limits defined at the critical points define 
distinct interacting field theories. For Nf = 6, however, simulations on a 16^ lattice 
with m = 0.01 give tentative evidence for a first order chiral transition [^, implying 
both that no continuum limit exists in this case, and that 4 < Nfc < 6. 

Both the nature of the transition and the value of Nfc are non-perturbative issues, 
inaccessible via either a standard perturbative expansion (which is non-renormalisable 
for d > 2), or a '\-/Nf expansion 0, ^. There have, however, been analytic attempts 
to investigate the transition via the truncated Schwinger- Dyson (SD) equations [^ 
H, 1^, 1^ . In the most systematic treatment , the SD equations are solved in ladder 
approximation in the the strong coupling limit; chiral symmetry breaking solutions 
are found for Nf < Nfc = 128/37r^ ~ 4.32. Note that in a continuum approach the 
value of Nfc manifests itself as the value of Nf below which a non-trivial solution 
can be found, eg. by bifurcation theory. In a later paper p|, Sugiura extended the 
solution to finite g^ for d G (2,4); for the case d = 3 his predictions for the critical 
coupling as a function of Nf read: 
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The essentially singular behaviour seen in (^) as Nf y Nfc is consistent with the 
existence of a conformal fixed point |p. Using a different sequence of truncations 
of the SD equations, however. Hong and Park ^ found chiral symmetry breaking 
solutions for all Nf, and predicted 
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which is manifestly non-perturbative in l/Nj. 

The differing predictions of the SD approach motivate the use of lattice field theory 
methods to address this problem. In this Letter we present results of simulations 
performed with values of Nf falling between the values Nf = 2,4 and 6 explored in 
previous works, in an effort both to map out the phase diagram in more detail, and 
to constrain further the value of Nf^. To start with, let us define the lattice action: 
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xi 9 xfi 
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where the indices i,j run over A^ flavors of staggered lattice fermion. The A^ are real 
auxiliary vector fields defined on the lattice links; gaussian integration over A^ yields 
a form of the action with explicit four-fermion couplings [Q, ^ (other variants of the 
lattice- regularised Thirring model are discussed in p|). Analysis of the spin- flavor 
content of staggered lattice fermions [l^ reveals that in 3 dimensions the number of 
continuum four-component physical flavors is given by Nf = 2N. Therefore use of 
the above action limits us to even Nf. If, however, the fermions are integrated out, 
then the resulting effective action for the A^ fields reads: 

Nf Nt 

Self = ^E4(^) - ^lndetA<[A]. (5) 

The effective action, though non-local, has an analytic dependence on A^j, and 
hence can be employed for odd or non-integer values.^ It can be simulated using a 
hybrid molecular dynamics algorithm, which evolves the {A} configuration through 
a fictitious time r by a combination of microcanonical and Langevin dynamics. In 
the limit of timestep Ar — > the {A} are distributed according to the equilibrium 



ensemble. In this work we have implemented the R algorithm of Gottlieb et al [11 

in this case the systematic errors are 0{N'^At'^) IjTT], 0. We found that on a 12^ 

lattice a value Ar = 0.01 was sufficient, with a mean interval r = 1.0 between 

refreshments. Our results typically arise from averages over 500 units of r. Fig. |I| 

shows a comparison between chiral condensate (xx) data obtained for Nf = 4, for 

^Note that the form of (0) permits even-odd partitioning, which means that we can simulate 
integer powers of detM, rather than detA^^A^, with a local action. 
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Figure 1: Chiral condensate (here denoted a) vs. m on a 12^ lattice for Nf = 4 with 

i/g' = 1.0. 



various values of the bare mass m, using both the hybrid algorithm and a hybrid 
Monte Carlo algorithm , which in principle is free from systematic error. 

In Sec. |2.1| we will present results from simulations with variable Nf at a fixed 
coupling l/g"^ = 1.0. The phase transition is located and found to fall at a value 
of Nf intermediate to the cases studied in earlier work p|, Q. Fits to a power law 
equation of state ansatz yield critical exponents consistent with this picture. In Sec. 
|2.2| we present results obtained with variable g"^ at a fixed Nf = 5, the goal being 
to determine if the transition remains continuous, or whether there is evidence for 
coexisting phases signalling a first order transition, as observed for Nf = 60]. Our 
conclusions, and the resulting phase diagram, are given in Sec. ^. 



2 Numerical Simulations 

2.1 Fixed 1/g^ = 1.0 

Using the algorithm as described above we performed simulations with variable Nf 
on a 12^ lattice at fixed 1/g'^ = 1.0, with bare mass m = 0.04, 0.03, 0.02. In Fig. || we 
plot the transverse and longitudinal susceptibilities Xt,i, using the definitions given in 
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Figure 2: Transverse and longitudinal susceptibilities vs. Nf on a 12'^ lattice with 
1/^2 = 10, m = 0.02. 

10, 1^. As described there, Xt can be calculated either from the integrated two-point 
pion propagator, or from (xx) via the axial Ward identity: the two methods give 
compatible results, the latter being the less noisy. These quantities are related to 
the inverse square masses of respectively the pseudoscalar and scalar bound states. 
For Nf > 3.6 they are roughly equal, indicating that the scalar and pion states are 
degenerate, and chiral sjTumetry is realised. For Nf < 2.8, in contrast, Xt ^ Xh 
indicating that the pion has become light as expected for a pseudo-Goldstone mode 
in a phase of broken chiral symmetry. Fig. ^thus offers evidence for a chiral symmetry 
breaking phase transition in the region Nf ~ 3.0. 

Fig. ^ shows results for the order parameter (xx) ^s a function of Nf, together 
with fits to a renormalisation-group inspired equation of state of the form 
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where 5, j3 are conventional critical exponents, and the parameter t expresses the 
distance from criticality. In previous work at fixed Nf, variable l/g"^ we have defined 
t = 1/ g^ — ^1 (}\\ here we define 



t = Nf- Nf,, 
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Figure 3: Chiral condensate a vs. A'^^ on a 12'^ lattice with 1/ g"^ = 1.0, showing fits 
to the equation of state @. The dashed line shows the fit in the chiral limit. 

where Nf^ denotes the value of Nf at the transition. Experience from previous studies 
|l|, ||, 1^ has shown that equations of state of the sort (|) provide adequate descriptions 
of the data in the vicinity of the transition, and that the resulting fitted exponents 
do not differ much from those extracted by more sophisticated finite volume scaling 
analyses. Also following previous work we fix (3 by the requirement 6 — 1/(3 = 1, in 
order to stabilise the fitting procedure; we refer to the resulting four parameter fit as 
fit I. The results obtained by fitting data for Nf G [3.0, 3.6] are given in Table |1|. Also 
shown is a five parameter fit II over the same range for which the constraint on j3 was 
relaxed. A number of comments are in order. First, the constraint of fit I appears to 
be well satisfied by the data, since from fit II the fitted value of 5 — 1/(3 = 0.97(8). 
Second, the fitted values of the exponents fall between those found for Nf = 2 and 
Nf = 4 from runs with variable l/g"^ [jl|, being numerically closer to the Nf = 4 case: 
the results for these from comparable fits are summarised in Table |^. This supports 
the idea that g1{Nf) is a critical line of fixed points along which critical exponents 
vary smoothly, with the numerical value of the exponents independent of the direction 
from within the ((7^, Nf) plane from which the critical line is approached. Third, as a 
note of caution, it should be noted that the criterion for the range of Nf over which 
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Table 1: Results for l/g'^ = 1.0 from fits to (|) on a 12^ lattice. 



Parameter 


Fit I 


Fit II 


Nf. 


3.060(30) 


3.060(44) 


A 


0.264(7) 


0.253(8) 


B 


3.9(9) 


4.4(1.3) 


5 


3.14(17) 


3.22(24) 


P 


— 


0.44(5) 


xVd.o.f. 


1.3 


1.0 



the fit was made was purely on the basis on minimising the x^ per degree of freedom, 
as in ||l], ^, 1^. Inspection of Fig. |^ suggests that the resulting fits may not describe 
the broken phase particularly well, and indeed fits which included smaller Nf values 
resulted in systematically larger values of 6. Our philosophy is to regard (^ merely as 
a phenomenological decription of the data; results from comparable fitting procedures 
applied to different models can thus be used to establish trends regardless of whether 
the true equation of state is ultimately of the form (|^) or not. 

Table 2: Comparison with critical parameters for the Nf = 2 and 4 models from a 

12'^ lattice from ref. [1]. 

Parameter Nf = 2 "iV/ = 3.06" Nf = A 

l/g'i (fit I) 1.94(4) 1.00(1) 0.66(1) 

6 (fit I) 2.68(16) 3.14(17) 3.43(19) 

P (fit II) 0.71(9) 0.44(5) 0.38(4) 



2.2 Fixed Nf = 5 

With the confidence that simulations with a non-integer number of staggered fiavors 
yield results consistent with our expectations for Nf ~ 3, we then conducted a series 
of simulations on 12^ lattices with similar parameters but this time with Nf = 5, a 
value which previous studies suggest is close to Nfc. Our results for (xx) are shown 
in Fig. |[ We have included results from simulations with m = 0.01, although it is 
likely that these are badly finite-volume affected; however, they do reveal evidence 
of discontinuous behaviour for 1/g^ ^ 0.5. Fits to the equation of state (|^), with 
the parameter t now given by 1/g^ — 1/ g1, were attempted, but were less successful 
than those of the previous section, in general either resulting in an unacceptably 
large x^ or not reaching convergence. In any case, the resulting fitted parameters are 
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Figure 4: Chiral condensate a as a function of 1/g^ for Nf = 5 on a 12^ lattice. 

strongly dependent on the fit range, and inclusion of the m = 0.01 points makes fit 
quality much worse. Our best results, coming from fitting data with m G [0.04, 0.02], 
l/g"^ G [0.4,0.6], are shown in Tab. |[ The fitted values of 6 and P are consistent 



Table 3: Results for Nf 


= 5 from fits to (^ on a 12'^ lattice. 


Parameter 


Fit I Fit II 



A 
B 
5 

xVd.o.i. 



0.437(4) 
1.15(3) 
190(40) 
5.6(1) 

2.7 



0.47(2) 

0.5(2) 

22(20) 

4.1(6) 

0.29(5) 

2.3 



with the trends of Tab. |^, although it is interesting to note that the fitted value of 
6—1/(3 = 0.6(2), indicating that in contrast to previous fits the data do not satisfy the 
usual constraint in this case. It could be argued that this disfavours the fit, since the 
combined assumptions of the degeneracy of scalar and pseudoscalar bound states in 
the chiral limit of the symmetric phase, plus the applicability of a power-law equation 
of state (I), naturally predict S — 1/ [3 =10, |13[- 
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Figure 5: Fisher plot of a^ vs. m/a for Nf = 5 on a 12^ lattice. 



To gain more insight we present the same data in the form of a Fisher plot (ie. 
(xx)^ vs. rn/ {xx)) ^^ Fig- H- This plot is devised to yield trajectories of constant 
1/g'^ which intersect the vertical axis in the chiral limit in the broken phase, and 
the horizontal axis in the symmetric phase. Departures from the mean field indices 
6 = 3, j3 = ^ are revealed as curvature in the lines, and departures from 6 — 1/(3 = 1 
are revealed by variations in the sign of the curvature between the two phases. If we 
ignore the m = 0.01 points, which are those closest to the horizontal axis, then the 
plot suggests a critical l/g^ ~ 0.5, with tentative evidence for 6 — 1/(3 < 1. Inclusion 
of the low mass points, however, suggests an accumulation of the constant coupling 
trajectories around a line which if continued would intercept the horizontal axis. This 
is similar to the Fisher plot for 12^ data for the model with Nf = 6, shown in fig. 6 of 
0, and suggests similarities between the two cases. For Nf = 6, tentative evidence 
for tunnelling between coexisting vacua on simulations with m = 0.01 on a 16^ lattice 
in the critical region was presented in [0], consistent with the chiral transition being 
first order. We tested this possibility by performing long Nf = 5 simulations on a 
16^ lattice for m = 0.01 and l/g"^ = 0.45, 0.47 (these runs requiring At = 0.002): a 
time history for the latter is shown in Fig. ^ Whilst the fluctuations are large, with 
excursions of 0(50%) about the mean, there is no convincing evidence for a two-state 
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Figure 6: Time history of a for Nf 



5, m 
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signal. Double gaussian fits to the histograms of the binned data, which revealed 
coexisting states via twin peaks in P|, were this time statistically indistinguishable 
from simple single gaussians. 

In summary, whilst there seems to be clear evidence for a chiral transition with 
1/(7^ = 0.50(5), the simulations performed to date are unable to determine the order. 
Fits based on the assumption of a continuous transition, whilst falling into the broad 
trend of existing data, are unsatisfactory and fail to satisfy the constraint 6 — 1/(3 = 1. 
A search for evidence of metastability consistent with a first order transition, on the 
other hand, proved negative. 

3 The Phase Diagram 

In this Letter we have used a fresh simulation algorithm to explore non-integer num- 
bers of staggered fermion flavors, and found results broadly consistent with earlier 
studies, namely that for Nf ~ 3 the transition is continuous with critical exponents 
intermediate between those of Nf = 2 and Nf = 4 models, and that for Nf = 5 
the order of the transition is difficult to determine, but shares features in common 
with the Nf = 6 transition, believed to be first order. This change in order of the 
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Figure 7: Phase diagram for the three dimensional Thirring model; the chirally broken 
phase is to the lower left. The dashed line shows a fit to @. 



transition with increasing Nj has also been observed in studies of QED4 |T^. The 
results also support the view that the partition function can be regarded as analytic 
in Nf, and that the non- locality of the fermionic action for Nf odd has no severe 
consequences. In this final section we collect together these new results with the best 
estimates obtained from previous finite volume scaling studies in [^, H, to plot the 
phase diagram of the Thirring model on the (l/f?^, Nf) plane in the chiral limit. The 
result is shown in Fig. |^. The shading of the points indicates the suspected order 
of the transition. The picture we have established is consistent with a critical line 
of fixed points, along which critical exponents vary smoothly, up to some Nfc — 5, 
whereupon the transition becomes first order. 

It is work comparing our findings to the analytic approximations of Refs. 
and p]. Certainly the curvature of our critical line 1/ g1{Nf) is consistent with both 
predictions @ and (^). A detailed comparison with the former, which predicts a con- 
formal fixed point at Nj = Nfc — 4.32 is hampered by the difficulties in identifying 
the strong coupling limit in the lattice model, since there is an additive renormalisa- 
tion relating I/qcont to l/g_ 



LATT ^^^ to non-conservation of the interaction current 
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in the lattice model [l|, 0. Therefore we are forced to identify Nf^ with that value 
at which the transition becomes first order, at which point no continuum limit exists 
for the lattice model. Our estimate for Nfc is thus roughly consistent; however, the 
trends in the exponents extracted from the lattice studies, summarised in Table. ^ 
do not support a conformal fixed point. At such a point we expect the exponent 5 to 
take the value 1, as for an asymptotically-free theory. The exponent r] is related to 



the anomalous dimension of the composite operator xx by 77 = d — 27^^ |T5[ , which 
in turn is predicted at the fixed point to be given by ||^: 

rf-2 



7xx 



2 

At a conformal fixed point r] should thus take the value 2. Finally, it is interesting to 

Table 4: Critical exponents as a function of Nf. 

Nf 6 7] 

2 2.75(9) 0.60(2) 

3.06(3) 3.14(17) — 

4 3.76(14) 0.26(4) 

5 5.6(1)? — 



note that a fit of the form 



7r2 



l/gl = A + Be^^\--Nsy (9) 

motivated by Eq. (|^) adjusted to allow for the renormalisation of 1/(7^, can be made 
to the data shown in Fig. |^; the fit is particularly successful if the Nf = 2 point 
is excluded, in which case the resulting coefficients are A = 0.28(2), B = 4.7(2), 
with x^/d.o.f.=0.6. This suggests that the scenario of [0, in which chiral symmetry 
breaking is predicted for all Nf and hence that Nfc = 00, cannot be excluded by the 
lattice studies to date. 
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